setwd("E:\\5hmc_file\\2_5hmc_yjp_bam\\ASM")
dir1="./bayes_pvalue_beta0/"
library(ggplot2)
library(dplyr)
group1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M50_M49","M28_M27","M30_M29","M26_M25","M35_M36","M18_M17","M20_M19","M22_M21","M40_M39")

for(i in 1:6){
  fn1=paste0(dir1,group1[i],".bayes_p.txt")
  file=read.table(fn1,head=T,sep = "\t")
  file3=file
  file3$FDR1=p.adjust(file3$normal_bayes_pvalue,method = "BH")
  file3$FDR2=p.adjust(file3$tumor_bayes_pvalue,method = "BH")
  file3$unitID=paste(file3$chrom,file3$position,file3$ref,file3$var,sep=":")
  bias=file3[file3$FDR1<0.1|file3$FDR2<0.1,]
  nobias=file3[file3$FDR1>0.1&file3$FDR2>0.1,]
  bias_normal=data.frame(VAF=bias$normal_var_freq,FDR=bias$FDR1,beta=bias$normal_bayes_beta0,lower=bias$normal_lower_beta0,upper=bias$normal_upper_beta0,unitID=bias$unitID)
  bias_tumor=data.frame(VAF=bias$tumor_var_freq,FDR=bias$FDR2,beta=bias$tumor_bayes_beta0,lower=bias$tumor_lower_beta0,upper=bias$tumor_upper_beta0,unitID=bias$unitID)
  bias_normal=bias_normal[bias_normal$FDR<0.1,]
  bias_tumor=bias_tumor[bias_tumor$FDR<0.1,]
  biasdata=rbind(bias_normal,bias_tumor)
  biasdata$VAF=as.numeric(gsub("%","",biasdata$VAF))/100
  
  nobias_normal=data.frame(VAF=nobias$normal_var_freq,FDR=nobias$FDR1,beta=nobias$normal_bayes_beta0,lower=nobias$normal_lower_beta0,upper=nobias$normal_upper_beta0,unitID=nobias$unitID)
  nobias_tumor=data.frame(VAF=nobias$tumor_var_freq,FDR=nobias$FDR2,beta=nobias$tumor_bayes_beta0,lower=nobias$tumor_lower_beta0,upper=nobias$tumor_upper_beta0,unitID=nobias$unitID)
  nobias_normal=nobias_normal[nobias_normal$FDR>0.1,]
  nobias_tumor=nobias_tumor[nobias_tumor$FDR>0.1,]
  nobiasdata=rbind(nobias_normal,nobias_tumor)
  nobiasdata$VAF=as.numeric(gsub("%","",nobiasdata$VAF))/100
  
  datatest=rbind(biasdata,nobiasdata)
  datatest$group="nobias"
  datatest[datatest$FDR<0.1,]$group="bias"
  p1=ggplot(datatest,aes(y=FDR,x=VAF,color=group))+geom_point(size=1)+
    scale_color_manual(values = alpha(c('red','#708090')))+theme_classic(base_size = 15)+geom_hline(yintercept = 0.1,color="red",linetype="dashed")+theme(legend.position = "none")
	fn=paste0("2021.01.09.bias.visualization/",group1[i],"_bias_visualization.pdf")
  ggsave(p,filename = fn,width = 5,height = 4)
}

#下面的代码用reads1和reads2画散点图，虽然能呈现偏移但是会暴露我们数据测序深度不够的短板
tmp1=data.frame(reads1=file3$normal_reads1,reads2=file3$normal_reads2,FDR=file3$FDR1)
tmp1$group=ifelse(tmp1$FDR<0.1,"significant","nosignificant")
tmp2=data.frame(reads1=file3$tumor_reads1,reads2=file3$tumor_reads2,FDR=file3$FDR2)
tmp2$group=ifelse(tmp2$FDR<0.1,"significant","nosignificant")
tmp=rbind(tmp1,tmp2)
ggplot(tmp,aes(y=reads1,x=reads2,color=group))+geom_point(size=1)+
  geom_smooth(data = subset(tmp,FDR>0.1),method = "lm",color="black")+
  geom_smooth(data = subset(tmp,FDR>0.1&reads1>reads2),method = "lm",color="red",se = TRUE)+
  geom_smooth(data = subset(tmp,FDR>0.1&reads1<reads2),method = "lm",color="red")+
  scale_color_manual(values = alpha(c('#708090','red')))+theme_light()

ggplot(tmp,aes(y=reads1,x=reads2,color=group))+geom_point(size=1)+
  scale_x_continuous(breaks=c(seq(0,75,25)),limits=c(0,75))+
  scale_color_manual(values = alpha(c('#708090','red',0.7)))+theme_classic(base_size = 15)

